(* ::Package:: *)

(* ::Input:: *)
Nw=5;
out=ConstantArray[0,{Nw,Nw}];
(*Func[x,y] = 1/Pi;*)
(*Func[x,y]=x / 1.9368054149716534;*)
(*Func[x,y]=(x^2-0.25) / 1.1097208122457478;*)
Func[x,y]=x*y/0.82239999999999935;
(*Func[x,y]=2x^2+2y^2-1;*)
(*Func[x,y]=3*x*x*x+3*x*y*y -2*x;*)
(*Func[x,y]=3*x*x*x+3*x*x*y -2*x;*)
acc=0
For[j=0,j<Nw,j++,
  jj=(j+.5)*2/(Nw)-1.0;
  y1=jj-1.0/Nw;
  y2=jj+1.0/Nw;
  For[k=0,k<Nw,k++,
    kk=(k+.5)*2/(Nw)-1.0;
    x1=kk-1.0/Nw;
    x2=kk+1.0/Nw;
    (*Print[{j,k,x1,x2,y1,y2}];*)
    result=NIntegrate[Boole[Sqrt[x*x+y*y]<1]*Func[x,y],
                   {x,x1,x2},{y,y1,y2},
                   PrecisionGoal->12];
    out[[1+j,1+k]]=result;
    acc+=result;
    Print[N[result,5]]
  ]
]

(*Print["Sum"]
Print[acc]
Print[out]
*)

Export["poly_xy.txt", out, "Table"]
